Для выполнения этого задания вам понадобятся данные о среднемесячных уровнях заработной платы в России:
В файле записаны данные о заработной плате за каждый месяц с января 1993 по август 2016. Если хотите, можете дописать в конец ряда данные за следующие месяцы, если они уже опубликованы; найти эти данные можно, например, здесь: http://sophist.hse.ru/exes/tables/WAG_M.htm
Необходимо проанализировать данные, подобрать для них оптимальную прогнозирующую модель в классе ARIMA и построить прогноз на каждый месяц на два года вперёд от конца данных.
In [1]:
#загружаем необходимые модули
%pylab inline
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
import warnings
from itertools import product
#обратное преобразование Бокса-Кокса
def invboxcox(y,lmbda):
if lmbda == 0:
return(np.exp(y))
else:
return(np.exp(np.log(lmbda*y+1)/lmbda))
In [2]:
#загружаем данные из файла
data = pd.read_csv('WAG_C_M_updated.csv', sep = ';', index_col=['month'], parse_dates=['month'], dayfirst=True)
data.rename(columns={'WAG_C_M':'wage'}, inplace=True)
print data.shape
data.head()
Out[2]:
Выведем данные на график:
In [3]:
plt.figure(figsize(15,7))
data.wage.plot()
plt.ylabel(u'Средняя зарплата')
pylab.show()
Проверка стационарности и STL-декомпозиция ряда:
In [4]:
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.wage).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.wage)[1])
По критерию Дики-Фуллера гипотеза о нестационарности ряда не отвергается.
Сразу видны следующие особенности данных:
Сделаем преобразование Бокса-Кокса для стабилизации дисперсии:
In [5]:
data['wage_box'], lmbda = stats.boxcox(data.wage)
plt.figure(figsize(15,7))
data.wage_box.plot()
plt.ylabel(u'Средняя зарплата после преобразования Бокса-Кокса')
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.wage_box)[1])
Размах дисперсии ощутимо уменьшился.
По критерию Дики-Фуллера нулевая гипотеза о нестационарности ряда не отвергается (0.72 > 0.05)
Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:
In [6]:
data['wage_box_diff'] = data.wage_box - data.wage_box.shift(12)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.wage_box_diff[12:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.wage_box_diff[12:])[1])
По критерию Дики-Фуллера можно отвергнеть гипотезу нестационарности (0.009 < 0.05), в остатках стало заметно меньше структуры, но полностью избавиться от тренда не удалось. Попробуем добавить ещё обычное дифференцирование:
In [7]:
data['wage_box_diff2'] = data.wage_box_diff - data.wage_box_diff.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.wage_box_diff2[13:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.wage_box_diff2[13:])[1])
Гипотеза нестационарности отвергается, остатки выглядят похоже на белый шум, и явного тренда больше нет.
Посмотрим на ACF и PACF полученного ряда:
In [8]:
plt.figure(figsize(15,12))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(data.wage_box_diff2[13:].values.squeeze(), lags=50, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(data.wage_box_diff2[13:].values.squeeze(), lags=50, ax=ax)
pylab.show()
Выбираем параметры нашей модели:
Q - значение последнего значимого сезонного лага на автокоррелограмме. Т.к. значимых лагов, кратных периоду (12) нет, то Q = 0
q - значение последнего значимого несезонного лага на автокоррелограмме. q = 1
P - значение последнего значимого сезонного лага на частичной автокоррелограмме. В данном случае это лаг = 12, поэтому возьмем значение P = 1
p - значение последнего значимого несезонного лага, меньшего величин периода, на частичной автокоррелограмме. p = 1
Начальные приближения: Q=0, q=1, P=1, p=1
In [9]:
#устанавливаем границы массивов наших параметров согласно начальным приближениям
ps = range(0, 2)
d=1
qs = range(0, 2)
Ps = range(0, 2)
D=1
Qs = range(0, 1)
In [10]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[10]:
In [11]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model=sm.tsa.statespace.SARIMAX(data.wage_box, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
Выводим удачные модели:
In [12]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
Лучшая модель:
In [13]:
print(best_model.summary())
Рассмотрим остатки модели:
In [14]:
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid[13:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])
Однако, уровни значимости критерия Льюинг-Бокса и критерия Стьюдента не слишком велики. Возможно, стоит перебрать больше значений параметров для поиска оптимальной модели?
In [15]:
ps = range(0, 4)
d=1
qs = range(0, 4)
Ps = range(0, 2)
D=1
Qs = range(0, 2)
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[15]:
In [16]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model=sm.tsa.statespace.SARIMAX(data.wage_box, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [17]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
Как видно, действительно удалось найти модель с лучшим рейтингом Акаики (24.27 < 36.75), хоть и более сложную
In [18]:
print(best_model.summary())
Рассмотрим остатки новой модели:
In [19]:
plt.figure(figsize(15,10))
plt.subplot(211)
best_model.resid[13:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])
Проверяем критерии:
Теперь посмотрим, как модель приближает данные.
In [20]:
data['model'] = invboxcox(best_model.fittedvalues, lmbda)
plt.figure(figsize(15,7))
data.wage.plot(label = 'фактические значения')
data.model[13:].plot(color='r', label = 'прогноз', linestyle = '--')
plt.ylabel(u'Средняя зарплата')
plt.legend()
pylab.show()
Визуально наша модель неплохо приближает реальные данные. Посмотрим прогноз.
In [21]:
#подбираем значения индексов для прогноза
best_model.predict(start=294, end=294 + 24)
Out[21]:
In [22]:
#Выводим прогноз
data2 = data[['wage']]
date_list = [datetime.datetime.strptime("2017-06-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,24)]
future = pd.DataFrame(index=date_list, columns= data2.columns)
data2 = pd.concat([data2, future])
data2['forecast'] = invboxcox(best_model.predict(start=293, end=294 + 24), lmbda)
plt.figure(figsize(15,7))
data.wage.plot(label = 'фактические значения')
data2.forecast.plot(color='r', label = 'прогноз', linestyle = '--')
plt.ylabel(u'Средняя зарплата')
plt.legend()
pylab.show()
Видно, что прогноз учитывает как сезонные колебания в данных, так и возрастающий тренд. Прогноз выглядит адекватно.